library(tidycensus)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggbeeswarm)
library(scales)
##
## Attaching package: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factor
library(ggiraph)
## Warning: package 'ggiraph' was built under R version 4.3.3
library(scales)
library(tigris)
## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.
us_income <- get_acs(
geography = "state",
variables = "B19013_001",
year = 2019,
survey = "acs1",
geometry = TRUE,
resolution = "20m"
)
## Getting data from the 2019 1-year ACS
## The 1-year ACS provides data for geographies with populations of 65,000 and greater.
## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
##
|
| | 0%
|
|====== | 9%
|
|============ | 18%
|
|================== | 26%
|
|========================= | 35%
|
|=============================== | 44%
|
|===================================== | 53%
|
|=========================================== | 61%
|
|======================================================= | 79%
|
|======================================================================| 100%
us_income_shifted <- us_income %>%
shift_geometry(position = "outside") %>%
mutate(tooltip = paste(NAME, dollar(estimate), sep = ": "))
gg_income <- ggplot(us_income_shifted, aes(fill = estimate)) +
geom_sf_interactive(aes(tooltip = tooltip, data_id = NAME),
size = 0.1) +
scale_fill_gradient(low = "lightgreen", high = "darkgreen", labels = label_dollar()) +
labs(
title = "Median household income by State, 2019",
caption = "Data source: 2019 1-year ACS, US Census Bureau",
fill = "ACS estimate"
) +
theme_void()
girafe(ggobj = gg_income) %>%
girafe_options(
opts_hover(css = "fill:cyan;"),
opts_zoom(max = 10)
)
Green was chosen for income to convey growth and prosperity, common associations with financial metrics. The legend shows dollar values for clarity, and the clean design, with minimal distractions, directs viewers’ attention to geographic patterns in income distribution across the states. Also, this interactive map could provide detailed value of each state when you click the state you want to check the income value.
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(ggplot2)
library(tigris)
options(repos = c(CRAN = "https://cloud.r-project.org"))
install.packages("ggspatial")
##
## The downloaded binary packages are in
## /var/folders/q7/ww3mcfxj1pv1181n2srlnm780000gn/T//RtmpCTNmiE/downloaded_packages
library(ggspatial)
library(ggrepel)
library(dplyr)
stations <- st_read("/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge/Highspeed_Stations/Highspeed_Stations.shp")
## Reading layer `Highspeed_Stations' from data source
## `/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge/Highspeed_Stations/Highspeed_Stations.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 74 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -75.35358 ymin: 39.90543 xmax: -75.07788 ymax: 40.11349
## Geodetic CRS: WGS 84
key_station_names <- c(
"30th Street Station", "69th Street Transportation Center", "Frankford Transportation Center",
"Norristown Transportation Center", "Fern Rock Transportation Center", "Arrott Transportation Center",
"NRG Station", "Chinatown", "8th", "Olney", "Girard", "Wynnewood Rd", "Ardmore Ave",
"Bryn Mawr", "Villanova"
)
key_stations <- stations %>%
filter(Station %in% key_station_names)
se_pa_counties <- counties(state = "PA", cb = TRUE, year = 2021) %>%
filter(NAME %in% c("Delaware", "Montgomery", "Philadelphia")) %>%
st_transform(st_crs(stations))
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|========================== | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 41%
|
|============================== | 43%
|
|=============================== | 44%
|
|======================================== | 58%
|
|================================================== | 72%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================= | 87%
|
|======================================================================| 100%
ggplot() +
geom_sf(data = se_pa_counties, fill = "grey90", color = "black", size = 0.5) +
geom_sf(data = stations, color = "blue", size = 1, alpha = 0.5) +
geom_sf(data = key_stations, color = "red", size = 3) +
geom_label_repel(data = key_stations, aes(geometry = geometry, label = Station),
stat = "sf_coordinates", nudge_y = 0.01, color = "darkred",
size = 3, min.segment.length = 0) +
labs(
title = "Key High-Speed and Major Stations in in Southeastern Pennsylvania",
caption = "Data source: SEPTA"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 18, face = "bold"),
plot.caption = element_text(size = 10)
)
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
## Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
This map highlights the key high-speed and major rail stations in Southeastern Pennsylvania, focusing on Philadelphia and its surrounding counties (Delaware, and Montgomery). The highlighted stations represent critical transit hubs and high-traffic stops, providing an overview of the region’s significant rail nodes. Key stations are marked with larger, red points to distinguish them from other stations, which are represented with smaller blue dots. Each highlighted station is labeled with a clear, bordered label, using ggrepel to prevent overlap and ensure readability.
hydro_features <- st_read("/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge/Hydrographic_Features_Arc/Hydrographic_Features_Arc.shp")
## Reading layer `Hydrographic_Features_Arc' from data source
## `/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge/Hydrographic_Features_Arc/Hydrographic_Features_Arc.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 8559 features and 28 fields (with 1 geometry empty)
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -75.46082 ymin: 39.80297 xmax: -74.89938 ymax: 40.24626
## Geodetic CRS: WGS 84
highlighted_features <- c("Schuylkill River", "Delaware River")
highlighted_hydro <- hydro_features %>%
filter(creek_name %in% highlighted_features)
highlighted_labels <- highlighted_hydro %>%
group_by(creek_name) %>%
summarize(geometry = st_union(geometry)) %>%
ungroup()
ggplot() +
geom_sf(data = se_pa_counties, fill = "grey95", color = "grey", size = 0.3, alpha = 0.5) +
geom_sf(data = hydro_features, color = "lightblue", size = 0.5) +
geom_sf(data = highlighted_hydro, color = "blue", size = 1.2) +
geom_label_repel(data = highlighted_labels, aes(geometry = geometry, label = creek_name),
stat = "sf_coordinates", nudge_y = 0.01, color = "darkblue",
size = 3.5, min.segment.length = 0) +
labs(
title = "Hydrographic Features in Philadelphia",
subtitle = "Highlighted rivers and waterways",
caption = "Data source: OpenDataPhilly"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 18, face = "bold"),
plot.subtitle = element_text(size = 14, face = "italic"),
plot.caption = element_text(size = 10)
)
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
This map visualizes the hydrographic network in Philadelphia, focusing specifically on the prominent waterways: the Delaware River and Schuylkill River. All hydrographic features in the city are displayed in a light blue to provide a contextual background, while the highlighted rivers are represented in a more vivid blue and labeled for clarity. By emphasizing the Delaware and Schuylkill Rivers, the map emphasizes the importance of these waterways to the city’s geography and infrastructure. This design allows viewers to quickly identify these key rivers, while still appreciating the complexity of Philadelphia’s overall hydrographic system.
planning_districts <- st_read("/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge//Vital_Natality_PD/Vital_Natality_PD.shp")
## Reading layer `Vital_Natality_PD' from data source
## `/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge/Vital_Natality_PD/Vital_Natality_PD.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2340 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -75.28031 ymin: 39.86747 xmax: -74.95575 ymax: 40.13793
## Geodetic CRS: WGS 84
planning_districts_filtered <- planning_districts %>%
filter(YEAR == 2023, METRIC_NAM == "count_births")
ggplot(planning_districts_filtered) +
geom_sf(aes(fill = METRIC_VAL), color = "white", size = 0.3) +
scale_fill_viridis_c(
option = "magma", direction = -1,
name = "Number of Births",
labels = comma
) +
annotate(
"text", x = -75.10, y = 39.85, label = "Lowest Births Here",
color = "red", size = 4, fontface = "bold"
) +
annotate(
"segment", x = -75.18, xend = -75.15, y = 39.9, yend = 39.85,
color = "red"
) +
labs(
title = "Philadelphia Vital Statistics - Natality (Births in 2023)",
subtitle = "Number of births by planning district",
caption = "Data source: OpenDataPhilly"
) +
theme_minimal() +
theme(
panel.grid = element_blank(),
panel.background = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 18, face = "bold"),
plot.subtitle = element_text(size = 14, face = "italic"),
plot.caption = element_text(size = 10),
legend.position = "right"
)
This choropleth map visualizes the number of births by planning district in Philadelphia for the year 2023. Each planning district is colored based on the total number of births, with darker shades indicating higher birth counts and lighter shades indicating lower birth counts. An annotation highlights the planning district with the lowest number of births. The goal of this map is to provide a clear and visually engaging way to understand natality trends across Philadelphia’s planning districts. By highlighting the lowest birth area, I want to form a narrative element, and inspire more public health, urban planning, or policy-making discussions.
install.packages("showtext")
##
## The downloaded binary packages are in
## /var/folders/q7/ww3mcfxj1pv1181n2srlnm780000gn/T//RtmpCTNmiE/downloaded_packages
library(showtext)
## Loading required package: sysfonts
## Loading required package: showtextdb
install.packages("ggtext")
##
## The downloaded binary packages are in
## /var/folders/q7/ww3mcfxj1pv1181n2srlnm780000gn/T//RtmpCTNmiE/downloaded_packages
library(ggtext)
landmarks <- st_read("/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge/Landmarks_AGOTrainingOnly-shp/8902d0d9-94eb-4763-a5e7-7dc4e412475a2020329-1-xj0ntw.uct4l.shp")
## Reading layer `8902d0d9-94eb-4763-a5e7-7dc4e412475a2020329-1-xj0ntw.uct4l' from data source `/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge/Landmarks_AGOTrainingOnly-shp/8902d0d9-94eb-4763-a5e7-7dc4e412475a2020329-1-xj0ntw.uct4l.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7740 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -75.28368 ymin: 39.85871 xmax: -74.95665 ymax: 40.13423
## Geodetic CRS: WGS 84
museums <- landmarks %>%
filter(FEAT_TYPE == "Museum")
showtext_auto()
font_add_google(name = "IM Fell English", family = "vintage")
museums_points <- museums %>%
st_centroid()
## Warning: st_centroid assumes attributes are constant over geometries
install.packages("ggimage")
##
## The downloaded binary packages are in
## /var/folders/q7/ww3mcfxj1pv1181n2srlnm780000gn/T//RtmpCTNmiE/downloaded_packages
library(ggimage)
museum_icon <- "/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge/MuseumIcon.png"
ggplot() +
geom_sf(data = landmarks, fill = "grey", color = NA, size = 0.1, alpha = 0.3) +
geom_image(
data = museums_points,
aes(x = st_coordinates(geometry)[, 1], y = st_coordinates(geometry)[, 2], image = museum_icon),
size = 0.02
) +
labs(
title = "<span style='font-size:16pt; font-family:vintage;'>Philadelphia Museums</span>",
subtitle = "<span style='font-size:12pt; font-family:vintage;'>Distribution of Museums in the City</span>",
caption = "Data source: OpenDataPhilly",
x = NULL,
y = NULL
) +
theme_minimal(base_family = "vintage") +
theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.background = element_rect(fill = "#F6F0E4", color = NA),
plot.background = element_rect(fill = "#F6F0E4", color = NA, size = 1),
panel.grid = element_blank(),
legend.position = "none",
plot.title = element_markdown(hjust = 0.5, margin = margin(b = 10)),
plot.subtitle = element_markdown(hjust = 0.5, margin = margin(b = 20)),
plot.caption = element_text(hjust = 0.5, face = "italic", size = 10, margin = margin(t = 20))
)
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
The IM Fell English font was chosen to evoke a classic, historical feel suitable for a map focusing on cultural landmarks like museums. A sepia-toned background (#F6F0E4) and muted black frame give the map a vintage and timeless appearance. The goal was to create a visually striking map that communicates the cultural richness of Philadelphia through its museums while adhering to the vintage theme.
crime_data <- read.csv("/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge/crime_incidents_philly.csv")
theft_data <- crime_data %>%
filter(!is.na(lat) & !is.na(lng) & text_general_code == "Thefts")
theft_sf <- st_as_sf(theft_data, coords = c("lng", "lat"), crs = 4326)
theft_sf <- st_transform(theft_sf, crs = 3857)
valid_x_min <- -75.3
valid_x_max <- -74.9
valid_y_min <- 39.8
valid_y_max <- 40.2
theft_sf_filtered <- theft_sf %>%
filter(
point_x != 0 & point_y != 0,
point_x >= valid_x_min & point_x <= valid_x_max,
point_y >= valid_y_min & point_y <= valid_y_max
)
ggplot() +
stat_bin_hex(
data = theft_sf_filtered,
aes(x = point_x, y = point_y),
bins = 30
) +
scale_fill_viridis_c(option = "inferno", name = "Theft Density", direction = -1) +
geom_label(
aes(x = -75.11, y = 39.95, label = "High Theft Area"),
fill = "white", color = "red", fontface = "bold", size = 4
) +
labs(
title = "Philadelphia Theft Incidents",
subtitle = "Density of theft incidents aggregated into hexagonal bins",
caption = "Data source: OpenDataPhilly"
) +
theme_minimal() +
theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 18, face = "bold"),
plot.subtitle = element_text(size = 14, face = "italic"),
plot.caption = element_text(size = 10)
)
This map visualizes the density of theft incidents in Philadelphia, by aggregating the data into hexagonal bins for spatial analysis. The color gradient, ranging from light yellow to deep purple, represents the density of theft incidents, where darker shades indicating higher densities. A prominent “High Theft Area” is highlighted with a red-bordered white label and an arrow pointing to the hexagonal bin with the highest density. Hexagons were chosen because of their aesthetic and analytical advantages, such as uniform adjacency and equal distance between centroids.For this map, I want to inform people the theft hot spot in Philadelphia and inspire public safety discussion.
library(sf)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
buildings <- st_read("/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge/LI_BUILDING_FOOTPRINTS/LI_BUILDING_FOOTPRINTS.shp") %>%
st_transform(crs = 4326)
## Reading layer `LI_BUILDING_FOOTPRINTS' from data source
## `/Users/luoxiaoyi/Desktop/2024 Fall/communication/30daysChallenge/30DayMapChallenge/LI_BUILDING_FOOTPRINTS/LI_BUILDING_FOOTPRINTS.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 544579 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -75.27943 ymin: 39.87237 xmax: -74.9573 ymax: 40.13778
## Geodetic CRS: WGS 84
buildings <- st_make_valid(buildings)
invalid_geometries <- sum(!st_is_valid(buildings))
if (invalid_geometries > 0) {
cat("There are still", invalid_geometries, "invalid geometries.\n")
} else {
cat("All geometries are valid!\n")
}
## All geometries are valid!
buildings <- buildings %>%
mutate(absolute_height = BASE_ELEVA + APPROX_HGT)
building_points <- st_centroid(buildings)
## Warning: st_centroid assumes attributes are constant over geometries
data <- st_coordinates(building_points) %>%
as.data.frame() %>%
mutate(height = buildings$absolute_height)
plot_ly(
data,
x = ~X, y = ~Y, z = ~height,
type = "scatter3d", mode = "markers",
marker = list(size = 3, color = ~height, colorscale = "Viridis", opacity = 0.8)
) %>%
layout(
title = "Philadelphia Building Heights",
scene = list(
xaxis = list(title = "Longitude"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "Height")
)
)
## Warning: Ignoring 6336 observations